#!/usr/bin/python
from numpy import *
from scipy import *
#from pylab import *
#Radiacion diaria

lat=40.416666667*pi/180
long=3.75*pi/180

#Entradas

#Radiacion Media de Madrid
Gdm=array([1.99 , 2.64, 4.32, 5.32, 6.28, 7.29, 7.47,  6.62 , 5.11, 3.4, 2.16, 1.72]) *1000

#Entradas

#Calculadas

Kt=Kd=Bd=ws=r=E0=Tau=dia=zeros((1,365), dtype='float')
# dias del ano

#meses:
Ene=linspace(1,31,31)
Feb=linspace(32,59,28)
Mar=linspace(60,90,31)
Abr=linspace(91,120,30)
May=linspace(121,151,31)
Jun=linspace(152,181,30)
Jul=linspace(182,212,31)
Ago=linspace(213,243,31)
Set=linspace(244,273,30)
Oct=linspace(274,304,31)
Nov=linspace(305,334,30)
Dic=linspace(335,365,31)

dia=linspace(1,365,365)

#declinacion solar (rad)
Tau=2*pi*(dia-1)/365 
decli=(0.006918-0.39912*cos(Tau)+0.070257*sin(Tau)-0.006758*cos(2*Tau)+0.000907*sin(2*Tau)-0.002697*cos(3*Tau)+0.00148*sin(3*Tau))

#Distancia Sol Tierra
E0=1.00011+0.034221*cos(Tau)+0.001280*sin(2*Tau)+0.000719*cos(2*Tau)+0.000077*sin(2*Tau)

#distancia entre el sol y la tierra en km
r=1.496E8/sqrt(E0)

#angulo de salida del sol
ws=-arccos(-tan(decli)*tan(lat))


#Radiacion Solar extraterrestre
Bd=24.0/pi*1367.0*E0*(cos(lat)*cos(decli))*(ws*cos(ws)-sin(ws))

Bdm_Ene=Bd[Ene.min():Ene.max()].mean()
Bdm_Feb=Bd[Feb.min():Feb.max()].mean()
Bdm_Mar=Bd[Mar.min():Mar.max()].mean()
Bdm_Abr=Bd[Abr.min():Abr.max()].mean()
Bdm_May=Bd[May.min():May.max()].mean()
Bdm_Jun=Bd[Jun.min():Jun.max()].mean()
Bdm_Jul=Bd[Jul.min():Jul.max()].mean()
Bdm_Ago=Bd[Ago.min():Ago.max()].mean()
Bdm_Set=Bd[Set.min():Set.max()].mean()
Bdm_Oct=Bd[Oct.min():Oct.max()].mean()
Bdm_Nov=Bd[Nov.min():Nov.max()].mean()
Bdm_Dic=Bd[Dic.min():Dic.max()].mean()

Bdm=([Bdm_Ene, Bdm_Feb,Bdm_Mar,Bdm_Abr,Bdm_May,Bdm_Jun,Bdm_Jul,Bdm_Ago,Bdm_Set,Bdm_Oct,Bdm_Nov,Bdm_Dic])


Ktdm= Gdm/Bdm
Kddm=([0,0,0,0,0,0,0,0,0,0,0,0])

for x in range(len(Bdm)):

	if Ktdm[x] <= 0.17 :
		 Kddm[x]=0.99
		
	elif Ktdm[x]<=	0.8 :
		Kddm[x]= 1.188-2.272*Ktdm[x]+9.473*Ktdm[x]**2-21.856*Ktdm[x]**3+14.648*Ktdm[x]**4
 
	elif Ktdm[x]>0.8 :
		Kddm[x]=0.225
	
#print Ktdm;
#print Kddm;


#print "dia=", dia[16]
#print "Tau=", Tau[16]
#print "decli=", decli[16]
#print "r=", r[16]
#print "ws=", ws[16]
#print "E0=", E0[16]
#print "Bd=", Bd[16]
print "Bdm=", Bdm
print "Ktdm=", Ktdm


#plot(Bdm)
#show()

Transicion_K=array([[0,	  1,	 2,	    3,	   4,	  5,	 6,		7,	    8,	  9	 ],
				[0.000, 0.300, 0.350, 0.400, 0.450, 0.500, 0.550, 0.600, 0.650, 0.700],
  				[0.300, 0.350, 0.400, 0.450, 0.500, 0.550 ,0.600, 0.650, 0.700, 1.000],
				[0.031, 0.058, 0.051, 0.052, 0.028, 0.053, 0.044, 0.085, 0.010, 0.319],
				[0.705, 0.694, 0.753, 0.753, 0.807, 0.856, 0.818, 0.846, 0.842, 0.865]])


K=Ktdm[11]

for x in Transicion_K[0,:]:
	  if Transicion_K[1,x]<=K<=Transicion_K[2,x]:
	  	Intervalo_K=linspace(Transicion_K[3,x],Transicion_K[4,x],11)

	  	
print Intervalo_K

if K<=1 :
	
		Mat_trans=array([[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000],
	 	 				[.167, .319, .194, .139, .097, .028, .042, .000, .014, .000], 	
	  					[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000],
	  					[.167, .319, .194, .139, .097, .028, .042, .000, .014, .000], 	
	  					[.250, .250, .091, .136, .091, .046, .046, .023, .068, .000], 
	  					[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000],
	  					[.167, .319, .194, .139, .097, .028, .042, .000, .014, .000], 	
	  					[.250, .250, .091, .136, .091, .046, .046, .023, .068, .000], 
	  					[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000]])

elif 2<=K<=3:
		Mat_trans=array([[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000],
	 	 				[.167, .319, .194, .139, .097, .028, .042, .000, .014, .000], 	
	  					[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000],
	  					[.167, .319, .194, .139, .097, .028, .042, .000, .014, .000], 	
	  					[.250, .250, .091, .136, .091, .046, .046, .023, .068, .000], 
	  					[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000],
	  					[.167, .319, .194, .139, .097, .028, .042, .000, .014, .000], 	
	  					[.250, .250, .091, .136, .091, .046, .046, .023, .068, .000], 
	  					[.229, .333, .208, .042, .083, .042, .042, .021, .000, .000]])


	#fila
for z in range(10):
		if Intervalo_K[z]<=K<=Intervalo_K[z]:
			fila=z

columna=0
Nrandom=rand()
	
#while K<Nrandom:
#		 columna=columna+1
#		k=K+ Mat_trans[fila,columna]
		
	#columna
		 
		 
		
